+(release notes: none)
+@amcastro-tri for feature review please. This PR is getting close to touch the core transfer algorithm, and I figure you might be interested in understanding how it works.

Reviewable status: LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

Working... I left a couple of questions.
Mostly I'd need some variables definitions in order to dig into the implementation of Sort().

Reviewed 5 of 8 files at r1, all commit messages.
Reviewable status: 8 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/spgrid.h line 28 at r1 (raw file):

/* A subset of SPGrid flags and configs we care about. */
struct SpGridFlags {

per this PR, this is struct and accessor below are not used.
Remove.

multibody/mpm/particle_sorter.h line 25 at r1 (raw file):

};

using Ranges = std::vector<Range>;

nit, I'd personally advocate for a more explicit name like say RangesVector to make sure it differentiates by more than one character from Range and also it comunicates quickly there's a std::vector underneath.
I'd also be okay with using std::vector<Range> all around.

multibody/mpm/particle_sorter.h line 32 at r1 (raw file):

 @pre data is sorted in ascending order and has at least two elements.
 @pre ranges != nullptr */
void ConvertToRanges(const std::vector<int>& data, Ranges* ranges);

can the input vector have repeated elements? say {1, 3, 3, 6, 9}? would that lead to an "empty" range? is that allowed?
Simply document requirements and/or add checks if possible.

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

 Given a set of MPM particles and the background grid, recall that we define the
 base node of a particle to be the grid node that is closest to the particle

Would this be correct?

Suggestion:

 Given a set of MPM particles and the background grid, recall that we define the
 base node of a particle to be the grid node that is closest (in the max norm) to the particle

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

  template <typename T, typename SpGrid>
  void Sort(const SpGrid& spgrid, double dx,
            const std::vector<Vector3<T>>& particle_positions) {

not sure if documented elsewhere, but could you say or reference what's the phisical space occupied by an SpGrid with grid size dx?
Do the grid's coordinates always go from 0 to kMaxGridSize*dx? i.e. the origin is always at the corner of the grid?

In other words, I guess the SpGrid defines a frame (definition?) and the particle_poisions are expressed in that frame?

multibody/mpm/particle_sorter.h line 216 at r1 (raw file):

  std::vector<int> sentinel_particles_;
  std::vector<int> data_indices_;
  std::vector<uint64_t> base_node_offsets_;

maybe I could guess from the impl, but I think a good idea to document here locally the intent of these individual variables. Is the comment at the top supposed to apply to sentinel_particles_ or to all of these three arrays?

Also,w hat's the difference between an "index" (as stored in data_indices_) and an "offset" (as stored by base_node_offsets_).
I know what an offset is from the docs of SpGrid, but I don't know what an "index" is.

Could you also provide a local definitiion for "sentinel particles"? or maybe reference the paper that provides that definition?

It'd help a lot with the review and I believe also with the future maintenance of this code.

multibody/mpm/particle_sorter.h line 220 at r1 (raw file):

   index of the particle. */
  std::vector<uint64_t> particle_sorters_;
  std::array<Ranges, 8> colored_ranges_;

ditto for these, probably local comments that explain their specific purpose would serve best those of us unfamiliar with this code, rather than a generic header for a group of variables.

multibody/mpm/grid_data.h line 116 at r1 (raw file):

  /* Sets all floating point data to zero, and set the index to be inactive.*/
  void set_zero() {

without reading the docs this methods sounds a lot like a harmless Eigen setZero().
IMO, the most "surprising" or really releveant effect of this function is to set to inactive, with the side effect of setting things to zero.

In other words, if the method is called set_inactive() I will probably be more prone to reading the docs rather than guessing this function's behavior (erroneously as just setting things to zero without changing index state).

Another option is to be very verbose and call this set_inactive_and_zero(), but I believe set_inactive() should suffice.

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/spgrid.h line 28 at r1 (raw file):

nvm

Reviewable status: 6 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

There are two reasons that this is named set_zero().

1. It's required by (and called from) the SpGrid class. See class documentation of SpGrid in spgrid.h.
2. It's important all floating point values are zero, because new values will later be added to this struct without resetting things to zero.

If you think the naming is confusing, I can be persuaded to rename this to reset and remove the existing reset function.

multibody/mpm/particle_sorter.h line 25 at r1 (raw file):

Done

multibody/mpm/particle_sorter.h line 32 at r1 (raw file):

Done

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

That would be correct, but I don't see the need to specify the norm. The result is the same with any reasonable norm (e.g. euclidean norm).

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):
SpGrid doesn't have the notion of dx. It's a lattice of integer valued 3 vectors.

Do the grid's coordinates always go from 0 to kMaxGridSize*dx? i.e. the origin is always at the corner of the grid?

This is abstracted away (as much as we could). As the user of SpGrid, you can basically treat it as an infinite lattice of Vector3i with a bijection from Vector3i to uint64t.

In other words, I guess the SpGrid defines a frame (definition?) and the particle_poisions are expressed in that frame?

I clarified the frame for particle positions.

multibody/mpm/particle_sorter.h line 216 at r1 (raw file):
All private member variables without documentation is explained in the accessor. I've made that clear.

Is the comment at the top supposed to apply to sentinel_particles_ or to all of these three arrays?

Just sentinel_particles.

Could you also provide a local definitiion for "sentinel particles"? or maybe reference the paper that provides that definition?

No paper discusses this concept AFAIK. This is an implementation detail. I made an attempt to clarify what a sentinel particle a bit more.

multibody/mpm/particle_sorter.h line 220 at r1 (raw file):

This is already explained in the accessor().

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

1. It's required by (and called from) the SpGrid class.

I know, but that's under our control, you can always change the function's name both here and in spgrid.h.

2. It's important all floating point values are zero

Yeah, I do see that also.

If you think the naming is confusing, I can be persuaded to rename this to reset

Well, you have a reset() now that sets things to NaN, I imagine for good reasons.

I think you could keep your reset() and have this method to have a more explicit funtion name. Something like zero_and_deactivate() could do?
I just think that within local scope (say in spgrid.h), saying set_zero() gloses over the important detail that now this data will be marked inactive.

You could also consider two separate calls if you hate the compound name, i.e. a set_zero() (sets all floating points to zero, leaves index untouched) and a set_inactive() (marks the index/data inactive).

multibody/mpm/particle_sorter.h line 38 at r1 (raw file):

ah, you are right. This is good. I extrapolated my mental picture of a box around the base node to a box associated with the max norm.

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

As the user of SpGrid, you can basically treat it as an infinite lattice of Vector3i

Ah, are you sayint that these i-j-k indexes can be negative also? and thus SpGrid mathematically spans all of R^3?

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

The order in which the particle data is preferrablly accessed.

Tryint to reconcilliate this with the doc for base_node_offsets().

Would it be true to say:

Returns a vector toaid the recursion of particle data in "sorted" order. This is the in which data is preferrably accessed when working with SpGrid.

I guess what I liked to get from this comment (mine might be wrong) is:

1. "Order", is this the "sorted order" referenced in base_node_offsets()?
2. "Prefered", is this a requirement for the best usage of SpGrid?

multibody/mpm/particle_sorter.h line 66 at r2 (raw file):

  template <typename T, typename SpGrid>
  void Sort(const SpGrid& spgrid, double dx,
            const std::vector<Vector3<T>>& particle_positions) {

Tryint to reconcilliate how indexes returned by data_indices() are assigned.
Are those indexes assigned based on the order in which particles are specified by particle_positions? i.e. the first position has index 0, the second index 1, etc.?

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

  std::vector<uint64_t> GetActiveBlockOffsets() const;

  /* The order in which the particle data is preferrablly accessed.

minor,

Suggestion:

preferablly

multibody/mpm/particle_sorter.h line 194 at r2 (raw file):

    for (int p = 0; p < num_particles; ++p) {
      const auto& particle_data = foo_data[data_indices[p]];
      // do something with particle_data.

ah, this is great. I guess what I didn't get from SpGrid's docs is that the offsets are not contiguous then? I was thinking I'd use offsets as 1D indexes. If true, maybe worth pointing out either here on SpGrid's docs.

multibody/mpm/particle_sorter.h line 209 at r2 (raw file):

   kernels, without write hazards between concurrent threads.

   Specifically, for each i from 0 to 7, colored_ranges()[i] is a collection of

I love this second paragraph, very clear.

Could you say something about this fixed (8) number of colors? why couldn't this be a different number? do we always have 8 colors and possibly empty ranges? or could we ever have more than 8 colors?
If needed, cite reference if you have it.

Reviewable status: 7 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

I could live without the current reset function and have the default state of the grid state to be zero and inactive. In that case reset is properly named.

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

Not technically R^3 because they are integer valued, but yes, they can be negative

multibody/mpm/particle_sorter.h line 66 at r2 (raw file):

I believe this resolved per our f2f, but feel free to continue the discussion here if you disagree.

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):
Working.

• "Order", is this the "sorted order" referenced in base_node_offsets()?

Yes

• "Prefered", is this a requirement for the best usage of SpGrid?

It's not neither a requirement nor necessarily absolutely the the "best", but it's the best I (and many before myself) can come up with.

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

Working

multibody/mpm/particle_sorter.h line 194 at r2 (raw file):

I think this is resolved from our f2f? Feel free to reopen and request changes if you disagree.

multibody/mpm/particle_sorter.h line 209 at r2 (raw file):

Working

Before we move on with the comments I was placing below, I have a higher level question. According to the comment below:

We sort particles first based on their base node offsets in the spgrid; if they are the same, then we sort the particles by their original indices.

That's a lexicographical order? It'd seem the bottleneck or most expesive part of Sort() is in std::sort.
I am wondering, wouldn't it suffice to create a std::vector<std::pair<uint64_t, int>> (a vector of {offset, index} pairs) and provide a custom coparison operator for std::sort()? In terms of implementation/documentaion it'd seem simpler, we no need to play with these custom masks along all the size checks to ensure correcntess. But I imagine somehow this would not do the trick for MPM applications?

Reviewable status: 10 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

oh, I see. I think that'd read nicely then.

multibody/mpm/particle_sorter.h line 59 at r1 (raw file):

right, thanks!

multibody/mpm/particle_sorter.h line 89 at r2 (raw file):

     because we cap the total number of allowed grid nodes. There are at most
     2^(3*log2_max_grid_size) number of grid nodes and that takes up
     3*log2_max_grid_size bits. The page bits and block bits have 64 - data bits

consider adding parentheses to ease the reading.
I had to scan this a couple of times untile I reliazed the hyphen was a minus sign.
Maybe use page_bits, block_bits and data_bits so that it looks more like math?

Ditto elsesewhere for consistency.

Suggestion:

have (64 - data_bits)

multibody/mpm/particle_sorter.h line 92 at r2 (raw file):

     in total, so the left most 64 - data bits - 3 * log2_max_grid_size bits are
     zero. So we left shift the base node offset by that amount and now we get
     the lowest 64 - 3 * log2_max_grid_size bits (which we name `index_bits`) to

Not the 64 -data bits - 3 * log2_max_grid_size? that is, what happened with the data_bits?

Code quote:

     zero. So we left shift the base node offset by that amount and now we get
     the lowest 64 - 3 * log2_max_grid_size bits (which we name `index_bits`) to

multibody/mpm/particle_sorter.h line 94 at r2 (raw file):

     the lowest 64 - 3 * log2_max_grid_size bits (which we name `index_bits`) to
     be zero. We use these bits to store the original index of the particles.
     With a typical log2_max_grid_size == 10, we have 44 bits to work with, more

I feel dumb. Shouldn't this be 34 - data_bits? (i.e. 64 - 3 * log2_max_grid_size = 34, minus the data bits)

Code quote:

we have 44 bits t

multibody/mpm/particle_sorter.h line 98 at r2 (raw file):

     bit unsigned integers which is enough to achieve the sorting objective. */
    const int index_bits = 64 - 3 * log2_max_grid_size;
    // TODO(xuchenhan-tri): This should be enforced at registration time.

Could you say more about what you mean with "registration" time?

Code quote:

 // TODO(xuchenhan-tri): This should be enforced at registration time.

multibody/mpm/particle_sorter.h line 111 at r2 (raw file):

        if constexpr (std::is_same_v<T, double>) {
          return particle_positions[p];
        } else if constexpr (std::is_same_v<T, float>) {

whey do we need to cast to double in the float case? is there a requirement I am missing?
That made me realize, shouldn't dx be of type T? if not, simply document why not.

multibody/mpm/particle_sorter.h line 121 at r2 (raw file):

          spgrid.CoordinateToOffset(base_node[0], base_node[1], base_node[2]);
      data_indices_[p] = p;
      /* Confirm the data bits of the base node offset are all zero. */

Could you say more about this assumption?
I'm confused about when you include the data_bits in the size computation and whether you don't.

Code quote:

 Confirm the data bits of the base node offset are all zero. 

Yes, that would be equivalent with the trick we are playing with masks.

However, as discussed in our f2f last week, the reason that we perform the masking trick is to eventually use non-comparative sort (e.g. radix) sort that has proven to be able to vastly reduce the time spent on sorting the particles compared to std::sort.

Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

Working

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

Done.

I added a note in the documentation of base_node_offsets() that clarifies base_node_offsets() and data_indices() uses the same particle order.

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

Done

multibody/mpm/particle_sorter.h line 209 at r2 (raw file):

Done.

I added a code example and a reference.

Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/grid_data.h line 116 at r1 (raw file):

Done.

TBH I don't recall this from our f2f, I think we only talked about the role of each indidual member in this class?
Anyway, that now makes sense. Thanks.

I'll coment below.

Reviewable status: 11 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

sorry, my correction was wrong also. Sitll misspelleed.

multibody/mpm/particle_sorter.h line 78 at r3 (raw file):

    /* We sort particles first based on their base node offsets in the spgrid;
     if they are the same, then we sort the particles by their original indices.

Per discussion thread above, please explain here right away why the need for this custom implementation of lexicographical ordering. (in retrospect the vector of pairs implementation could've been a good first step for this PR, but I won't block you on that).
That is, explain right away your intetion to use a radix sort.

Code quote:

    /* We sort particles first based on their base node offsets in the spgrid;
     if they are the same, then we sort the particles by their original indices.

multibody/mpm/particle_sorter.h line 133 at r3 (raw file):

    // TODO(xuchenhan-tri): use a more efficient sorting algorithm. This can be
    // a bottleneck for large number of particles.
    std::sort(particle_sorters_.begin(), particle_sorters_.end());

upate this TODO to mention the intention to use radix sort, which essentially justifies the complex piece of code above.

Reviewable status: 11 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

multibody/mpm/particle_sorter.h line 89 at r2 (raw file):

Done

multibody/mpm/particle_sorter.h line 92 at r2 (raw file):

Recall that the bit representation of the base node offsets look like:

page_bits | block_bits | data_bits

The data bits are on the right. So it has nothing to do with the amount of left shifting we will perform.

I also reorganized the explanation of how the bit manipulation game is played. Hopefully it helps.

multibody/mpm/particle_sorter.h line 94 at r2 (raw file):

Yes, it should be 34 (64 - 30 = 34, not 44). But it's 34, period, not (34 - data_bits). This is because data_bits are zero too! We are free to overwrite it with useful information.

multibody/mpm/particle_sorter.h line 98 at r2 (raw file):

Done.

multibody/mpm/particle_sorter.h line 111 at r2 (raw file):

Here we are doing some computation where we want to stripe autodiff to get the underlying floating point scalar type.

Ideally, we want to keep float and double and convert autodiff to double, but I haven't figured out a clean way to do that, so instead, I converted everything to double.

multibody/mpm/particle_sorter.h line 121 at r2 (raw file):

This assumption is guaranteed by SPGrid. The 64 bit representation of all grid nodes (i.e. the "offset") have zero data bits.

multibody/mpm/particle_sorter.h line 185 at r2 (raw file):

Done

multibody/mpm/particle_sorter.h line 78 at r3 (raw file):

Done. I moved up the TODO below.

multibody/mpm/particle_sorter.h line 133 at r3 (raw file):

Done, I moved this TODO up a bit.